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Abstract — A generalized state space representation of a dynamical 
system with random modes is presented. The new formulation includes 
a term, in the dynamics equation, which depends on the most recent 
state's Unear minimum mean squared error (LMMSE) estimate. This can 
be used to model the behavior of a feedback control system featuring 
a state estimator. The measurement equation is allowed to depend on 
the previous LMMSE estimate of the state, which can be used to 
represent the fact that measurements are obtained from a vaUdation 
window centered at the predicted measurement and not from the entire 
surveillance region. The matrices comprising the system's mode constitute 
an independent stochastic process. The proposed formulation generalizes 
several problems considered in the past, and allows a unified modeling 
of new ones. The LMMSE optimal filter is derived for the considered 
problem. The approach is demonstrated in the context of target tracking 
in clutter and is shown to yield performance that is competitive to that 
of several popular nonlinear methods. 

Index Terms — State estimation, target tracking, fault detection and 
isolation, clutter and data association, hybrid systems. 



I. INTRODUCTION 

State estimation in dynamical systems with randomly switching 
coefficients is an important problem in a variety of applications. 
The most natural examples ai^e maneuvering target tracking and 
fault detection and isolation (FDI) algorithms, featured, e.g., in 
aerospace navigation systems. The standard modeling presumes that 
the dynamics of the continuously-valued state, and, possibly, its 
measurement equation, are controlled by an evolving mode that takes 
discrete values. This is the well known concept of hybrid systems fTl . 

Various problems have been formulated using a hybrid state 
space framework. In problems involving uncertain, or intermittent 
observations, such as l2|-j5J, the mode affects the matrices of the 
measurement equation. In maneuvering target tracking applications, 
considered in, e.g., (6|-(9l, the mode usually affects the matrices of 
the dynamics equation. 

We consider a state space representation of dynamical systems 
with randomly switching coefficients, that is more general than the 
commonly used one in two major aspects. First, we allow the system 
input to depend on the latest estimate of the state, as is common 
practice in closed loop control systems. In this work, the state 
estimate is taken to be the linear minimum mean squared error 
(LMMSE) estimate. Likewise, we allow the measurement equation to 
depend on the latest state estimate (taken as well to be the LMMSE 
estimate). This can be used to represent the fact that observations 
are not taken in the entire admissible space, but rather in a small 
validation window set about the predicted measurement of the state. 

It is well known (7) that, even for the simplest case of indepen- 
dently switching modes, the optimal estimate of the state cannot 
be obtained without resorting to exhaustive enumeration. Therefore, 
significant efforts have been dedicated to developing suboptimal 
approaches for state estimation in hybrid systems and especially 
for the important subclass of jump linear systems (JLS). Origi- 
nally derived in (TJ, the generalized pseudo-Bayesian (GPB) filter 
approximates the conditional expectation of the state of a Markov 
JLS (MILS) by pruning the measurement history. Perhaps the most 
successful suboptimal nonlinear method is the interacting multiple 
model (IMM) algorithm, proposed in (8j, in which history pruning 
is accompanied by a wise cooperation of the primitive Kalman 
filters comprising the algorithm, thus improving the performance 



over GPB with only a slight increase in computational compexity. 
Alternatively, one may consider optimality within a narrower family 
of linear filters. While (2|, (3) and others considered special cases 
of problems involving random modes, De Koning IIOI derived a 
Kalman filter (KF) like algorithm for the most general case of JLS 
with independent random modes. Costa llll derived a linear optimal 
scheme for the most general case of a standard MILS by means 
of state augmentation. Recently, it was shown in |12| that, in some 
cases, parts of the state may be estimated optimally while others in 
a linear optimal manner. 

In this paper we consider optimality within the family of linear 
filters for a generalized state space representation of dynamical 
systems. One of the attractive properties of LMMSE estimators is 
their minimax optimality 1131 . That is, they attain the lowest worst- 
case MSE, in comparison to any other (even nonlinear) method, 
for the same first- and second-order moments of the underlying 
distributions. We derive a linear optimal filter that may be con- 
veniently implemented in a recursive form, eliminating the need 
for unbounded memory, and without state augmentation. Our filter 
reduces to previously reported results when the parameters of the 
underlying problem are appropriately adjusted. 

As an illustration of the approach, we show how the problem of 
tracking a target in clutter may be formulated within our framework. 
We show that in this setting, the performance of our filter is 
competitive to that of several classical nonlinear methods. 

The remainder of this paper is organized as follows. In Section HH 
we define and discuss the proposed state space modeling approach 
and survey some related contributions. The linear optimal recursive 
state estimation algorithm is developed in Section [in] An application 
of the approach to target tracking in clutter, followed by a representa- 
tive numerical study, is presented in Section HVl Concluding remarks 
are provided in Section [Vl 

II. System Model and Related Work 



We consider the dynamical system 



Xk+i = AkXk + BkUk + CkWk 
yk = HkXk + GkVk + FkXk-i, 



(la) 
(lb) 



in which Xk £ R" and yk G are the state and measurement 
vectors at time k, respectively. The noise processes {wk} and {vk} 
constitute zero-mean unity-covariance white sequences, and xo is a 
random vector with mean xq and second-order moment Po- 

We consider two versions for the modeling of the input vector 
itfc. In the first (standard) case, Uk is a known deterministic input. 
However, as in some cases Uk serves as a closed loop control signal, 
it is common practice to let it depend on the most recent estimate 
of the state. Thus, in the second variant we set Uk = Xk, where 
Xk is the LMMSE estimate of Xk using the measurement history 

yk = {yi,-- ■,yk}- 

Likewise, the term Xk-i in the measurement equation is the 
LMMSE estimate of Xk-i based on the measurement history Vk-i- 
Affecting the measurement at time k, the term FkXk-i can be 
used to represent the fact that observations are not taken in the 
entire (feasible) space, but, rather, in a small (admissible) validation 
window, set about the predicted measurement. 

Finally, the system mode, Mk = {Ak, Bk, Ck, Hk,Gk, Fk}, 
constitutes an independent random process whose distribution at time 
k, p{A4k), is known. The random quantities {wk}, {vk}, {Mk}, 
and Xo are assumed to be mutually independent. 

We seek to obtain the linear optimal estimate Xk+i using the 
measurements 3^fc+i. It will be shown in the sequel that Xk+i in 
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our setting conveniently possesses thie recursive form 

Xfc+i = LfcXfc + KkVk+i + JkUk- (2) 

Namely, the desired linear optimal solution will be shown to fuse 
the (extrapolated) previous estimate with the most recently acquired 
measurement vector, effectively eliminating the need to store the 
entire measurement sequence. For the case in which Uk — Xk, the 
terms LkXk and JkXk in l|2j may be grouped together. 

Note that the described problem does not require the system 
mode to assume values in a discrete domain. In addition, the above 
formulation allows evolution not only of the entries of the matrices 
constituting the mode, but also of their dimensions. This observation 
provides a convenient basis for treating problems that, to the best of 
the our knowledge, have not been previously considered, as will be 
discussed in detail in Section [TVl 

For the setting in which no feedback terms are present (namely, 
Bk = and Fk = 0), several variants and special cases of the 
presented problem have been considered in the past. Nahi devised 
a linear optimal estimator having the recursive form ^ for the prob- 
lem of state estimation in the presence of multiplicative measurement 
faults. In this setting, only the matrix Hk in ^ is random. A gen- 
eralization of Nahi's work has been recently proposed in ||5j, where 
multiplicative faults were accompanied by additive ones representing 
measurement biases. This can be cast within the formulation l[2j, 
where only the matrices Hk and Gk are random. De Koning 1101 
considered the most general case of independently switching modes 
and Costa II II developed, by means of state augmentation, a recursive 
implementation of the LMMSE filter for systems with modes obeying 
Markov dynamics taking values in some discrete domain. Several 
additional contributions related to the considered problem include (3), 
in which a generalized version of Nahi's problem was addressed 
by allowing correlated fault indicators; (4), that allowed correlations 
between subsequent fault variables; j6j, that proposed a linear optimal 
estimator for the static multiple model (SMM) problem II14I ; and 1151 , 
that considered linear optimal estimation in systems with bounded 
and unbounded numbers of actuator faults. Nonlinear suboptimal 
solutions for related problems were proposed in (7)-(9), 1161 and 
references therein. 

III. Linear Optimal Recursive Estimation 

We now derive the LMMSE recursive filter for the case where Uk 
is a known deterministic signal. The result for the stochastic case is 
given in Section UlI-EI 

Let Yk be the random vector obtained by concatenating the 
measurements comprising . The desired recursive form ^ will be 
derived using the following lemma, which follows directly from ll71 
p. 190] and from the linearity of the MMSE estimator in the Gaussian 
case. 

Lemma 1. Let x, y and z be random vectors and let x{z) and 
x{y,z) denote, respectively, the LMMSE estimates of x using z, and 
using both y and z. In addition, let y{z) denote the LMMSE estimate 
of y using z. Then 



x{y,z) = x{z) +r^yry^y, 



(3) 



where y ~ y ~ vi^) <^'''d ^xy is the cross-covariance matrix between 
the random vectors x and y. 

Letting z = Yk and y = yk+i and using Lemma [T] we obtain the 
following form for the LMMSE estimate of Xk+i using Vk+i'- 



where Xf.^-^ is the LMMSE estimate of Xk+i using the measurements 
yk = {yi, yk}, yk+i = yk+i - y^+i^ ™d K+i the LMMSE 
estimate of yk+i using yk- It follows that 

f fe+i = E [Ak] £fc + E [Bk] Uk (5) 
i/fe+i = E [Hk+i] ifc+i + E [Fk+i] Xk 

= {¥.[Hk+i]^[Ak]+^[Fk+i])xk+^[Hk+i]^[Bk]uk. (6) 

Plugging (O in we identify the desired matrix coefficients Kk, 
Lk, and Jk of ^ as follows: 



(7) 



Lk = {I- KkE [Hk+i])E [Ak] - A-fcE [Fk+x] (8) 
3k = {I-KkE[Hk+i])E[Bk\. (9) 

We now compute the covariance terms ^xf^^iy^.^^ and '^y^^^iin^^i- 

A. Computation of Tx^^^y^^i 
Since is unbiased, 

r^fc+iSfc+i = E [(sfc+i - E [xk+i]){yk+i - Vk+iV] 

^E[xk+i{yk+i~yk+^V], (10) 
Substituting dlbt and 



^Xk+i{Hk+iXk+i + Gk+iVk+i + -Ffc+i£fc)^j 
- E [sfe+i((E [Hk+i] E [Ak] + E [Fk+i])xkf] 
-E[xk+i{^[Hk+i]E[Bk]ukf] . (11) 



r - — F. 



Utilizing the independence of Xk+i and Vk+\, and canceling out 
identical terms i ll It becomes 

r^,+i«.+i = E[xk+ixl+{\E[H^+^] - E[xk+ixl]E[Al]E[Hj+i] 
-E[xk+i]ulE[Bl]E[Hl+,]. (12) 

Before proceeding, we define the following matrices 

T.k=E[xkxl] (13) 

Afc ^ E[xkxl] = E[xkxJ] (14) 

Tk=E[xk]ul =E[xk]ul (15) 

Afe = UkuJ , (16) 

where the RHS of l ll4t follows from the orthogonality principle 
satisfied by Xk, the estimation error of Xk, and the RHS of l |15| l from 
the unbiasedness of Xk. Note that E^, A^, and are symmetric. 
Now, utilizing the independence of Xk and Wk, 

E ^Xk+ixl^ = E |^(Afea;fc + BkUk + CkWk)xX^ 

= E \{AkXk + BkUk)xl^ 

= E[Ak]Kk^'E[Bk]Tj, (17) 

which yields for l ll2t 

^^k+ivk+i = Sfe+iIEi-fffc+i] 

- {E[Ak]Ak +E[Bk]T]:)E[Al]E[H^+,] 
-E[xk+i]ujE[Bj]E[H^+,]. (18) 

To express E [a;fc+i] in terms of E [xk] we substitute i fTat 



Xk+l — ifc+l +^X^^ly^,+ lTy|^^-^y|^^-^yk + l, 



(4) 



E [xk+i] = E [AkXk + BkUk + CkWk] 
^E[Ak]E[xk]+E[Bk] Uk. 



(19) 
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Using the latter in l llSt and rearranging the summands, the final 
expression for rx^_f_iy^.^i reads 

^-.+iy,+i = ^k+iml+i] - (E[ylfel(AfeE[^^] + TfcE[i3j]) 

+ E[Bfe](TlE[Aj] + AfcE[Bj]))E[//J+i]. (20) 



B. Computation of Vy^^-^y^^-^^ 

Recall that y'^jf_-^ is the LMMSE estimate of yu+i and, conse- 
quently, j/fe+i is orthogonal to y^J^j. Hence, 

^Vk+ivk+i =E[(2/fc+i -j}fe+i)yj+i] 

= ^[yk+iyl+i] - ^[Vk+iVk+i] 

= E[yfc+iyJ+i] - (E[//fe+i]E[A-] +E[Ffe+i])E[ifej/J+i] 
-E[//fe+i]E[Bfc]^fcE[yJ+i], (21) 

where we have utilized the independence of {xk,yk+i} and 
{j4fe, -Bfe, //fc+i, Jfe+i}. Using l llbt and the independence of 

{ife,a;fe+i}, {_fffe+i,Gfc+i,i^fe+i} and v^+i, we have that 

E[ifej/J+i] = E |^ifc(i/fc+ia;fc+i + Gk+iVk+i + 
= E \xk{Hk+iXk+i + Ffc+iife)^j 
= E[xfca;;r+i]E[//J+i] + AfeE[FfeT^i], (22) 
which, using l |17t , becomes 

E[ifcyJ+i] = Afe(E[^^]E[i/J+i] +E[F,7+i]) 

+ TfcE[Bj]E[//,T+i]. (23) 

Due to the independence of {xk+i,S:k}, ffe+i, and {Hk+i, Gk+i} 

^[yk+iyl+i] 

= ^[{Hk+iXk+i + Gfc+iUfc+i + -Ffc+iife) 

X [Hk+ixk+i + Gfc+iUfc+i 
= ¥.[Hk+iXk+ixZ+^Hj+-^] + E[Gfe+iUfe+ii;J+iGfc+i] 
+ ^[Fk+iXkxl F^+^] + E[//fe+ia;fc+ii;Ji?'jr+i] 
+ E[Ffe+ia;fea:J+i7y,T+i]. (24) 

Consider, for example, the term E \j^k+iXkx'[j^_iH^j^_^. From the 
smoothing property of the conditional expectation, 

¥.[Fk+iXkxJ+iHk+i] = E \E[Fk+iXkxl+iHk+i \ Fk+i,Hk+i\ 

= E |^Ffc+iE[ifca;fc+i | Fk+i, Hk+i\Hj+i 

= E \Fk+Mxkxl+,\H^+^] , (25) 

where the conditioning on Hk+i and Fk+i was omitted due to the 
independence of {xk+i,Xk} and {Hk+i,Fk+i}. 

Similarly, bearing in mind that E [a::fc4.ia;^_^_j] = Efe+i, 
E [vk+ivj^i'j = J, and E [ifeX^T] = A^, we obtain: 

nHk+iXk+ixJ+,H^+i] = W.[Hk+iT.k+iH^+i] (26) 
E[Gfe+i«fc+iwJ+iGj+il = E[Gfc+iGj+i] (27) 

(28) 



E[Fk+iXkxJ fJ+,] = E[Ffc+iAfeFjr+i]. 

For future reference, we also note that 

E[AkXkxlAj] = E[AkJ:kAj] 
E[AkXkul bJ] = ElAkTkBj] 
E[BkUkujBj] = E[BkAkBj] 
E[CkWkwJcJ]=E[CkCj]. 



(29) 
(30) 
(31) 
(32) 



Substituting ^j} in l|25ll, and using llU-lllU in l|24ll 

E[yk+iyl+i] = E[Hk+i^k+iHj+i]+E[Gk+iG^+i] 
+ E[Fk+iAkF^+j] 

+ E [Hk+i{^ [Ak] Ak+E [Bk] )F,7+i] 

+ E [Fk+i{^kE[Aj] + TfcE[Bj])iJ,T+i] . (33) 

In addition, we obtain in a straightforward manner 

E [yk+i] = (E [Hk+i] E [Ak] + E [Fk+i])E [xk] 

+ E[Hk+i]E[Bk]uk. (34) 

Using ll26ll, 1(27), and l|28j in and substituting l|23ll, and OU 
in J2U we obtain the final expression for ^y^^iyk+i- 

^vk+ivk+i = ^[Hk+i'^k+iH^^i] +E[Gfe+iGfc+i] 
+ E[ffe+iAfeffeVi] - E[ffe+i]AfeE[ffeVi] 

- E [Hk+i] ElAk]{AkE[Aj] + rkE[Bj])E[Hj+,] 

- E [Hk+i] E[BkKrjE[Aj] + AkE[Bj])E[Hl+,] 
+ E [i7;fc+i(E[Afc]Afe + E[Bk]r]:)F^+,]^ 

+ E [Fk+ii^kE[Aj] + TkE[Bj])Hj+,] 

- E [Hk+i] E [Ak] AkE[F^+-,] - E [Fk+i] AfcE[Aj]E[//J+i] 

- E[Ffc+i]TfcE[Bj]E[iJ,T+J-E[//fc+i]E[Sfc]T^E[F,7+i]. (35) 

C. Computation of the Second-Order Moments 

The second-order moment matrix of Xk+i is given by 

Efe+i = E[xk+ixJ^i] 

= E ^{AkXk+BkUk +GkWk){AkXk+BkUk +CkWk)^^ 



- E[ylfeSfe Aj] + E[AfcTfcBj] + E[BfcTj 
+ E[BkAkBj]+E[CkCi^], 



(36) 



where we have utilized the independence of Xk, Wk and 
{Ak,Bk,Gk}, and 
Next consider A^+i. Using ^ we have: 

Afc+i = E |^(LfcXfe + Kkyk+i + JkUk)xJ+i^ 
= (Lfc + KkE[Fk+i])E[xkxJ+i] 
+ KkE [Hk+i] Efc+i + JkUkE[x 

fc+i]- 

Using l |17) the latter becomes 

Afc+i = Lk{AkE[Aj] + TfcE[Bj]) + Jk{rjE[Aj] + AkE[Bj]) 
+ Kk{EFk+i{AkEAj + TkEBj) + E[Hk+i]Ek+i). (37) 

Finally, 

Tfc+i = E[a:fe+i] uj+i. (38) 
Note that is known and does not need to be computed recursively. 

D. Algorithm Summary 

The complete algorithm comprises the following steps: 

Initialization: xq = xo. So = Po + xqXq , Ao = xqXq , To — 

XqUq , Ao = UouJ . 



SUBMITTED TO IEEE TRANS. ON AUTOMATIC CONTROL 



4 



Algorithm 1 



Input: yk+i, itfc+i, x^, E[xk], E^, A^, T^, 
1: Compute E [Ak], E [Bk], E [CkC^] 

E [AkTkBj] , and E [Bfc A^BJ] . 
2: Compute E [sfc+i] and E^+i using Eqs. il9\ and l |36t 



E[Ffc+i], 
and 



Tafc+iSfc+i using Eqs. (|20ll and OS. 



3: Compute ^E[i?fc+i], E [Gfe+iGfc+i 

E [Hfe+i(E [Ak] Afc + E [Bfc] Tj)Ffer_^i 
4: Compute F, 

5: Compute Jffe, L^, and using Eqs. Q, ([8](, and ([9l(, 
6: Compute A^+i and T^+i using Eq. ( 1371 ) and dSSI l. 
7: Compute £fc+i using Eq. (O. 
Output: Xk+i, E[xk+i], Efe+i, Afc+i, T^+i 



Recursion: For fc = 1, 2, . . . perform the routine of Algoritiim[T] 
Notice that since the distribution of the mode at each time is known, 
the expectations of steps [T] and [3] of the algorithm may be calculated 
in a straightforward manner. For example, if the mode variables take 
discrete values, these expectations reduce to probability-weighted 
summations of the corresponding matrices. However, in some cases, 
as demonstrated in Section ITVl closed form expressions exist for the 
above expectations. 

E. Random Inputs 

In the second variant of JTat . in which Uk = XkAt turns out that the 
roles played by Ak and Bk are identical. Specifically, after replacing 
life with Xk, at each step of the derivation of SectionUn] Ak and Bk 
are multiplied by the same quantities. Thus, the filter for the modified 
problem is obtained from the one described in Alg.[T]by replacing Ak 
with Ak + Bk and nullifying Uk and Tfe. An alternative derivation, 
based on the orthogonality principle, may be found in 1181 . 

IV. APPLICATION TO Target Tracking in Clutter 

In this section we demonstrate the proposed concept by casting the 
classical problem of tracking in clutter within our formulation, and 
applying the LMMSE filter of Section HiH 

A. System and Clutter Models 

Consider a single target obeying a linear dynamical model. The 
evolution of the state is obtained from i fTat by setting Ak — A, 
Bk ~ 0, and Ck ~ C, resulting in 



Xk+l 



Axk + Cwk 



(39) 



Here A and C are deterministic matrices, accounting for the state 
dynamics and process noise covariance, respectively. 

At time k, the target state is observed via the the linear equation 



true 

Vk 



™ I true 



(40) 



where v]^"^" represents measurement noise. In addition to the mea- 
surement y^k^°, at each time, a number of clutter detections are 
obtained. These will be denoted as {y^'i}^^^, where N is total 
number of detections. Clutter measurements originate from false (or 
ghost) targets and do not carry any information about the target of 
interest. They are, however, indistinguishable from true detections. At 
each time, the clutter measurements are assumed to be independent of 
each other, of the clutter measurements at other times, and of the true 
state and observation. In addition, we assume that these measurements 
are uniformly distributed in space. As an illustration, we present in 
Fig. [Da single realization of a one dimensional scenario in which the 




400 600 

Time 



1000 



Fig. 1: True target position (solid line) and all the obtained measure- 
ments (pluses) vs. time. 



true target position is accompanied by the actual measurements and 
clutter observations. 

To correctly model the distribution of the clutter detections, we 
note that typically, at each scan, the sensor initiates a validation 
window centered at the predicted target position, and the algorithm 
processes only those measurements obtained within the window. 
Since the clutter detections are uniformly distributed in space, they 
are also uniformly distributed within the validation window. 

We define the measurement vector yk at time k to be the concate- 
nation of all measurements. A'' — 1 of which correspond to clutter, 
and one originating from the true target. The location of the true 
measurement within this concatenated vector is, of course, unknown 
to the algorithm. This setting can be modeled using \lhi by letting 
the mode A^^ be distributed as 



Mk = {Hk,Gk,Fk} 



(41) 
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where diag (^Ai A2 ■ ■ ■ Am) stands for a block diagonal ma- 
trix. Here, Gei is the square-root of the covariance matrix associated 
with the clutter, which, as mentioned above, is uniformly distributed 
in the validation window. 

The first realization of {Hk, Gk, Fk} in ( I41t . for example, corre- 
sponds to the scenario in which the first of the A'^ observations is the 
true target measurement, yjj^'', generated according to i40i . while the 
other A'^ — 1 measurements are clutter, each of which is generated 
according to 



el 
VkA 



= HnomAXk-1 + GclUfc.i, i = 2, . 



(42) 



Here, HnomAxk^i is the predicted true measurement at time k, 
which is also the center of the validation window set by the sensor, so 
that clutter measurements acquired at time k are uniformly distributed 
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around this quantity. The overall number of measurements in the 
validation window, A^, is assumed to be known, but may vary in 
time. Thus, the dimensions of Hk, Gk, and may depend on k. 

Notice that we assumed, for simplicity, that the true measurement is 
always present in the validation window. To account for the possibility 
that the true measurement does not fall in the validation window, the 
option 



{Hk,Gk,Fk} = {0,In GciAn <» HnomA} , 



(43) 



needs to be added to set of possible realizations of {Hk, Gk, Fk} 
appearing in | |4U . Here, the symbol ® stands for the Kronecker 
product, Ijv is an A'' X 1 vector comprising all ones, and In is 
the N X N identity matrix. 

Finally, notice that if A'^ = 0, namely, there are no measurements 
in the validation window, then ^ becomes, at the absence of Uk, 
Xk+i ~ LkXk, which corresponds to a simple prediction (time up- 
date) without consecutive measurement update, as could be expected. 

B. Matrix Computations 

To invoke the algorithm presented in Section [in] we need to 
compute the application-dependent terms outlined in Steps [T] and [3] 
of Alg. [T] Although these may be evaluated numerically, via direct 
summations, in the present example closed-form expressions exist, 
as we show next. 

Since the matrices of the dynamics equation are deterministic, 
E [Ak] =A,E [Bk] = 0, E [CkCn = GG^, E [AkTkB^] = 0, 
E [BkAkBj] = 0, and E [Ak'SkAl] = AEkA'^ . Further, accord- 
ing to the probability distribution defined in ( 14 It . 



E[Hk + l] = -^IjV ® Hnom 

N ~ 1 

E [Fk+i] = ^ In ® Hnon,A. 
Computation of the remaining expectations yields 



(44a) 
(44b) 



E[Hk + lT.k + lHj^l] — — /jV (g) //nomSfc + l//Joni, (45) 

E[Gk+iGj+^] = ^/iv ® (CnomGL^ + (iV - 1)GciGJ) , (46) 
and 

where S is defined by 



(47) 



N = 1 



^{{N^2)1n1J, + In), N>1. 



(48) 



Finally, 

E 



[Hk+i{E [Ak] Ak + E [Bk] rJ)F^+,'\ = E [Hk+iAAkF^+,'\ 

nom / • 

(49) 



= -^(livljv - In) ® (^HnomAAkA iJ„ 

Following our assumption, the spatial distribution of the clutter 
measurements is uniform in the validation window. Thus, GciGJi is 
the covariance matrix of a random vector uniformly distributed in 
the validation window. For computational simplicity, we assume that 
the window is rectangular, so that the covariance computation boils 
down to finding the variance of scalar uniform random variables. 

Note that Bk = implies Jk = 0, and together with GciGj?j, 
the above moments allow the computation of the linear optimal filter 
coefficients, Kk and Lk, according to l|7j and l[8j, respectively. 



C. Discussion 

It is easy to see that, in the present case, F 



Vk+lVk+l 



is a block 



diagonal matrix such that the blocks along its main diagonal are 



D = - 

J 

4 

In addition, 

^^k+lVk+l 



Hnon-iAAkA _ffn 



+ -J^ Hnon-iT,k + lH2 



^G 



N - 1 
N 



GciGc 



1 



AAkA 



N 



{T.k+1-AAkA' 



■■■HJ 



Therefore, the filter gain, Kk, in ^ becomes 

J^k — L aifc + iafc+i-'- Vk + iVk + l 



N 



(Jlk+i- AAkA^)Hj,^D- 



(50) 



(51) 



(52) 



Bearing in mind that yk+i is a concatenated vector of all the 
observations acquired at time k + 1, the product Kk.yk+i in (O 
is nothing but the average of the measurements constituting yk+i, 
pre-multiplied by (E^+i — AAkA^)Hj^^-^D^^ . Consequently, the 
linear optimal estimator for tracking a target in clutter is a KF-like 
algorithm, operating on the average of all detections in the validation 
window. In this respect, its mode of operation resembles that of other 
classical methods. For example, the probabilistic data association 
(PDA) 1191 method applies a KF on a weighted average of all 
measurements in the window. The weights are inversely proportional 
to the exponent of the Mahalanobis distance between each of the 
detections and the predicted measurement, corresponding to the 
covariance of the innovation. Similarly, the nearest neighbor (NN) 
technique 1201 , 1211 applies a KF on the measurement nearest to the 
predicted estimate. This operation can be thought of as a weighted 
average in which one detection is assigned a weight of 1, while the 
rest are assigned weights. 

D. Numerical Study 

In this section we demonstrate the performance of the linear 
optimal filter for tracking a target in clutter by comparing it with that 
of the NN and PDA methods. We consider a one-dimensional tracking 
scenario, in which a target is represented via a two dimensional 
state comprising position and velocity information, Xk ~ {pk, Vk)^ ■ 
Starting at xq = (0, 0)^ with Po — (o o)' the target is simulated 
for 1000 time units by running the dynamical equation i39\ with 



The true measurement is generated by computing j40l t with 



f^nom — (1 0), Gnon 



0.32. 



(54) 



For a (one-dimensional) window of length d, the clutter variance 
of ^ is GciGj = d'^/12. 

We compare the performance of our filter with that of the NN and 
PDA filters. All the considered algorithms are designed to reduce the 
mean-squared estimation error in either a heuristic or an analytical 
manner. However, when dealing with tracking in clutter, using the 
MSE as the only performance measure may result in misleading 
conclusions. This is due to the fact that once the estimated state draws 
away from the true measurement, clutter measurements are likely to 
be treated as the true ones, eventually resulting in target loss. When 
this happens, the algorithms' MSE becomes meaninglessly large. 
Therefore, the MSE should be treated as a meaningful performance 
measure only as long as the target is not lost. We use two measures 
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Fig. 2: Position RMSE vs. clutter density. 
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Fig. 3: Track loss times vs. clutter density. 



of performance to evaluate the performance of the algorithms. The 
first is the time until the target is lost, defined as the third time the 
distance between the predicted position and the true state has deviated 
by more than five standard deviations of the (true) measurement noise. 
The second measure is the average squared error calculated over the 
time interval until the first of the three algorithms has lost track. This 
makes the comparison fair in the sense that none of the algorithms 
incorporates meaninglessly large errors corresponding to a lost target. 
A good algorithm is expected to have long track loss times, while 
maintaining low average squared errors. 

We test the algorithms versus an increasing clutter density. To this 
end we define p to be the average number of clutter measurements 
falling in an interval of one standard deviation of the (true) measure- 
ment noise. Averaged over 1000 independent Monte Carlo runs, the 
average squared position enws and track loss times, versus p, are 
plotted in Figures |2l and[3l respectively. 

It is readily seen that the linear optimal filter attains the longest 
track loss time while keeping the estimation enws at a reasonable 
compromise between the nonlinear PDA and NN algorithms. 



the state. These additional terms may be used to represent closed- 
loop control input, and measurement validation window, respectively. 
We derived a linear optimal recursive algorithm for this setting, and 
illustrated the proposed approach in the context of target tracking in 
clutter. In this situation, the new filter demonstrates competitive per- 
formance, when compared with several classical nonlinear methods. 
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V. CONCLUSION 

We proposed a general formulation of dynamical systems with 
independently switching coefficients, where the dynamics and mea- 
surement equations are allowed to depend on previous estimates of 



